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Abstract 

A numerical method based on Matrix Product Formalism is proposed to study the 
phase transitions and shock formation in the Asymmetric Simple Exclusion Process 
with open boundaries and parallel dynamics. By working in a canonical ensemble, 
where the total number of the particles is being fixed, we find that the model has 
a rather non-trivial phase diagram consisting of three different phases which are 
separated by second-order phase transition. Shocks may evolve in the system for 
special values of the reaction parameters. 
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During the last decade the study of one- dimensional driven lattice gases has 
been of increasing interest because besides their application in different areas 
of non-equilibrium physics they show a variety of fascinating properties such 
as non-equilibrium phase transitions and spontaneous symmetry breaking [1]. 
They have also let us study the one-dimensional shocks i.e. discontinuities in 
the density of particles on the lattice over a microscopic region. Depending on 
the dynamics these models can be divided into two different classes: models 
with sequential (continuous time evolution) and parallel (discrete time evo- 
lution) dynamics. These models can also have open boundaries (where the 
particles can enter or leave the lattice) or periodic boundary conditions (with 
conservation of the number of particle). 

Different approaches have been used to study the one-dimensional out-of- 
equilibrium models. Among these approaches the Matrix Product Formalism 
(MPF) is the one which allows us to study these models under sequential 
or parallel dynamics with different boundary conditions [2]. According to the 
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MPF the stationary probability distribution function of the system can be 
written in terms of the matrix element (for the open boundary case) or trace 
(for the periodic boundary case) of products of non-commuting operators. 
Several models have been studied using the MPF. In some cases they can 
be solved exactly; however, in most of the cases exact solutions can only be 
found in small regions of the reaction parameters space. In a recent work [3] 
it has been shown that the MPF can be used in order to numerical study of 
one-dimensional reaction- diffusion models with sequential dynamics and con- 
servation of the number of particles. In the present letter we aim to show that 
this approach can also be applied to the numerical study of one-dimensional 
reaction-diffusion models with parallel dynamics and conservation of particles. 
The Asymmetric Simple Exclusion Process (ASEP) with parallel dynamics will 
be considered as an example. 

In the ASEP with parallel dynamics ^ classical particles move on a one- 
dimensional lattice of length L (which is assumed to be an even number) with 

open boundaries. The bulk dynamic is deterministic and consists of two half 
time steps. In the first half time step particles at even sites move one step to 
the right provided that their rightmost sites are empty. In this step both the 
first and the last sites are also updated. From the first (last) site the particles 
are injected (extracted) with the probability a {/3) provided that the target 
site is empty (occupied). In the second half time step only the odd sites are 
updated and the particles at these sites move to the right in the same way. 
The parallel updating scheme is shown in Fig. 1. The ASEP with open bound- 
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Fig. 1. Parallel updating scheme for the ASEP with open boundary. 

aries and parallel dynamics has been originally proposed in [4]. Later it was 
studied using the MPF in [5] . They have both worked in the grand canonical 
ensemble where the total number of particles in the system is not conserved. 

The result is that in the large-L limit the model has two different phases: a 
high-density phase for ct > /3 and a low-density phase for a < A first-order 
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phase transition also takes place at the transition point a = /3 where the den- 
sity profile of particles changes linearly along the lattice. The linear density 
profile of particles is related to the superposition of shocks in the system. In 
order to study the shocks one may consider the model on a ring (a lattice 
with periodic boundary conditions) either in the presence of a slow particle 
[6,7,8], which is sometimes called an impurity, or a slow link [9]. Equivalently 
one can leave the boundaries open and restrict the total number of particles 
by working in a canonical ensemble. In the present letter we adopt the lat- 
ter scenario i.e. we leave the boundaries open; however, the total number of 
particles on the lattice M (and therefore their density p — ^) is restricted to 
being a constant. 

It is shown that the stationary probability distribution function of the ASEP 
with parallel dynamics can be written as [5] 

L 

P{n, ...,rL)^^{W\ l[[d,,_M\V) (1) 

in which we have defined 

Oi := TiD + (1 - Ti)E , Oi := nD + (1 - Ti)E (2) 



and Ti G {0, 1} is the occupation number at site i. The denominator in (1) is 
a normalization factor. The matrices (D, E) and (D, E) are square matrices 
and besides the vectors \V) and (yV\ are acting in an auxiliary space which 
might have either a finite or an infinite dimension. They also stand for the 
presence of particles and holes at odd and even sites respectively and should 
satisfy the following quadratic algebra 

[E, E] = [D,D] = , ED = [E, D] , DE ^ 

{W\E{1 -a) = {W\E , {W\{aE + D) = {W\D (3) 
(1 - P)D\V) = D\V) , {E + PD)\V) = E\V). 



It has also been shown that (3) has a two-dimensional representation for a ^ (3 
[5] 
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The normalization factor Z^, in (1), which will be called the grand canonical 
partition function for a > (3 hereafter, can easily be calculated using the fact 
that J2{c} ■ ■ ■ ) tl) = 1. It is found 

Zl = {W\[{E + b){E + D)]i\V) = (1 - - (1 - (5) 

For a < (5 the grand canonical partition function should be defined as —Zl- 
Let us investigate the phase transitions of the model in this case using the 
classical Yang-Lee theory [10,11]. One can apply the Yang-Lee theory by find- 
ing the roots of the partition function of the system as a function of one of 
its intensive variables. It can readily be seen that in the thermodynamic limit 
L — > oo the zeros of (5) as a function of a lie on a circle of radius (3 (a). 
This predicts a first-order phase transition at a = /3. Since the representation 
of the algebra (3) has a finite dimensional representation, we expect that the 
density profile of the particles has an exponential behavior in both phases 
a < /3 and a > (3. This is in agreement with the results obtained in [4,5] for 
the phase diagram of the model. As we mentioned above at the phase transi- 
tion point one finds a linear profile for the density of particles on the lattice. 
This is related to the superposition of shocks which can be anywhere on the 
lattice. In what follows we will investigate the phase diagram of the model 
with fixed number of particles by introducing a canonical partition function 
as 

L k 

Zl,M- E 5{M -Y,n){W\X{[^2^-l02^]\V) (6) 
{ri=0,l} 1=1 i=l 

in which M and L arc the number of particles and the length of the system 
respectively and 5{x) is the ordinary Kroncckcr delta function (5^ q- The oper- 
ators Oi and Oi are also defined in (2). It can easily be verified that (6) can 
be written as 

Zl,m = {W\Coef ficient[{CC)^ ,M]\V) (7) 

in which we have defined 
C ■=E + x]D 

(8) 

C:^E + xD 

where x is a free parameter and Coef ficient[Expr, n] gives the coefficient of 
in the expression Expr. Note that since Expr is a matrix then Coef ficient[Expr, n] 
gives the coefficient a;" of each of its entries. A similar formula has already 
been introduced in [3] which can be used for the numerical calculation of the 
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canonical partition function of one-dimensional reaction-diffusion models with 
sequential dynamics. The formula (7) can be useful for numerical studying of 
phase transitions in out-of-equilibrium systems with parallel dynamics. For 
the ASEP with parallel dynamics we will calculate (7) as a function of a and 
f3 for finite values of L and M by using the representation of its quadratic 
algebra (4). 

Now we investigate the density profile of the particles on the lattice which is 
defined as 

T,{c}^[n, - ■ ■ 



where C is any configuration with a fixed number of particles and P(rj, ■ ■ ■ , r^) 
is given by (1). The formula (9) can be written in terms of the operators 
D, E, D, E and the vectors \V) and however, since the operators for the 
existence of particles at even and odd sites are not the same we find different 
expressions for the density of particles at odd and even sites. It can easily be 
seen that the density of particles at odd and even sites can be written as 

{W\Coef ficient[{CCy-\ k]b Coef ficient[C{CC)^-\ M - k - 1\\V) 



and 



{W\Coefficient[{CCy-^C, k]DCoef ficient[{CC)^-\ M - k - l]\V) 



respectively where i = l,---,^. The formulas (7), (10) and (11) are quite 
general and can be used for any one-dimensional reaction-diffusion model with 
open boundaries, parallel dynamics and only one class of particles. 
In order to study the phase transitions of the ASEP in canonical ensemble 
with open boundaries and parallel dynamics, we use the classical Yang-Lee 
theory by investigating the zeros of (7) in the complex plane of both a and P 
for different values of L and M. The particle concentration at odd and even 
sites will also be calculated from (10) and (11) in each phase. Our numerical 
calculations show that in the thermodynamic limit {L — > oo, M — > oo,p ~ 
^) the model has two different phase diagrams depending on the density of 
particles on the lattice p: For p < \ the phase diagram of the model consists of 
a low-density and a shock phase which are separated by a second-order phase 
transition at /3c = 2p. The reason that the phase transition is of a second-order 
is that the zeros of (7) as a function of /3, approach the positive real-/3 axis 
at an angle ^ [12]. The canonical partition function (7) as a function of a, 
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does not have any real and positive root smaller than one in this case. In the 
low-density phase f3 > Pc the density of particles at both even and odd sites 
changes exponentially along the lattice; however, in the shock phase P < /3c as 
we move from the first to the last site of the lattice the particle concentration 
changes abruptly from /? to 1 for even sites and from to 1— a for odd sites. For 
p > I the phase diagram of the model consists of a high-density and a shock 
phase which are separated by a second-order phase transition at Oc = 2(1 — p). 
As for the case p < | the zeros of (7) as a function of a, approach the positive 
real- a axis at an angle | and therefore the phase transition is of second-order 
[12]. The canonical partition function (7) as a function of f3 does not have 
any real and positive root smaller than one in this case. In the high-density 
phase a > etc the density of particles at both even and odd sites changes 
exponentially; however, the shock phase in this case takes place at a < etc 
and the structure of the particles concentrations in quite similar to the shock 
phase for p < |. The phase diagrams for both p < | and p > | cases are 
shown in Fig. 2. At p = | the phase diagram of the model consists of only 
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Fig. 2. The phase diagrams of the ASEP in canonical ensemble with parallel dy- 
namics and open boundaries for p < ^ (left) and p > \ (right). The small curves 
show the behaviors of the density profile of particles in each phase. 

a shock phase independent of the values of a and (3. For fixed values of the 
injection and extraction probabilities (cto, /3o) the phase diagram of the model 
consists of three phases which are determined by the density of particles p. For 
< p < Y we are in the low-density phase in which the density of particles 
changes exponentially along the chain. For ^<p<l — ^we are in the shock 
phase. Finally, for 1 — ^<p<lwe are in the high-density phase in which 
the density profile of particles is an exponential function. 
The phase diagram of the ASEP in canonical ensemble with parallel dynamics 
and open boundaries can also be obtained by studying the grand canonical 
partition function of this model defined as 

Zl{x)^{W\{CC)^\V) (12) 
in which C and C are given by (8) and x plays the role of the fugacity of 
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particles. This can easily be calculated and we find 



Zl{x) = (1 - /3)Q;^+ixt(/? + x{l - - (1 - - a) + xa)^ 

+a(l-a)/3(l-/3)(l-x)(a:Q; + /3) (13) 



Now let us study the total density of particles as a function of the fugacity x 
p{x)^j-^\ogZ,{x). (14) 



In Fig. 3 we have plotted p{x) as a function of a; for L = 15000, a = 0.8 and 
f3 — 0.4. As can be seen, the density of particles is an increasing function of 
X up to a critical point Xq where a finite discontinuity takes place. Above the 
critical point xq the density increases very slowly and remains smaller than 
unity unless a — ^ 0. Our numerical calculations show that the discontinuity 
starts from pi = f and ends at p2 = 1 — f (see Fig. 3). This is quite in 
agreement with the picture that we got for the phase diagram of this model. 
For X < Xo and x > xq the density of particles is determined by their fugacity; 

p(x) 0.l[ ^^^^^^^^^^^^^^^^^^^^^^^^^^ . 
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Fig. 3. The density of particles p{x) as a function of their fugacity x for L = 15000, 
a = 0.8 and /3 = 0.4. The critical densities are pi = f = 0.2 and p2 = 1 — f = 0.6. 

however, at xq the fugacity does not fix the density. Here is where we have 
shocks in the system. This phenomenon has already been observed in other 
models where the Bose-Einstein condensation takes place and the conservation 
of the number of particles is broken [13]. 

In this letter we have investigated the phase transitions and shock formation in 
the ASEP with open boundaries and parallel dynamics in a canonical ensemble 
where the total number of particles is kept fixed. We have found that the phase 
diagram of the model depends of the value of the density of particles on the 
system. The system can be in any of its three accessible phases: the low-density 
phase, the high-density phase and the shock-phase. In the shock-phase, the 
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shock position is fixed and determined by the number of particles while in the 
hight-dcnsity and the low-density phases the density profile of the particles at 
both even and odd sites have exponential behaviors. Since the representation 
of the associated algebra (3) is finite dimensional, we do not expect the density 
profile of particles on the lattice to have algebraic behavior [14]. 
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